###############################################################
## The following code is part of the replication archive of  ## 
## Genovese, Kern & Martin (2016) "Policy Alternation", ISQ  ## 
###############################################################

#install.packages("systemfit")
#install.packages("stargazer")
#install.packages("nls2")
#install.packages("visreg")
library(foreign)
library(maps)
library(mapdata)
library(nlme)
library(nls2)
library(lattice)
library(systemfit) 
library(foreign)
library(stargazer)
library(coefplot)
library(effects)
library(colorspace)
library(splines)
library(ggplot2)
library(visreg)
library(systemfit) 
#library(reshape2)
#library(readstata13)

##################
# Figure 2: Maps #
##################

map('worldHires', c('Germany', 'Austria', 'Czech','Poland'), xlim=c(5,25), ylim=c(45,55))
text(10, 51, "Germany", cex=1.55)
text(14.5, 47.5, "Austria", cex=1.55)
text(19.5, 52.5, "Poland", cex=1.55)
text(16.5, 49.5, "Czech R", cex=1.55)

###################################
# Close Countries: Read the data  #
###################################

# Read the data:
scr <- read.dta("~/Dropbox/.../Replication/green taxes and subsidies/OECD subsidies lag.dta")
scr$ccode<-as.factor(as.character(scr$ccode))
scr$year<-as.factor(as.character(scr$year))

##########################################################
# Close countries, OLS w/ unstandardized rows: Table A.2 #
##########################################################

tax.raw1 <- lm(greentax_pc_cons  ~  lag_gtax_pc_cons  + l_subsidy_rgdppc   + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp  +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax + ccode  , data=scr)
tax.raw2 <- lm(greentax_pc_cons  ~  lag_gtax_pc_cons  + l_subsidy_rgdppc   + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp   +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  + l_sl_neigh_subsidy_rgdppc_raw + ccode   , data=scr)
tax.raw3 <- lm(greentax_pc_cons  ~  lag_gtax_pc_cons  + l_subsidy_rgdppc   + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp   +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax + l_sl_neigh_subsidy_rgdppc_raw + I( l_sl_neigh_subsidy_rgdppc_raw*actual_economic_flows) + ccode    , data=scr)
summary(tax.raw1)
summary(tax.raw2)
summary(tax.raw3)

###############################################
# Close countries, 3SLS w/ GDP p.c.: Table A.4 #
###############################################

tax2   <- greentax_pc_cons  ~ lag_gtax_pc_cons +    l_subsidy_rgdppc   + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party  + l_sl_neigh_subsidy_rgdppc +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax #   + ccode 
subsidy2  <- subsidy_rgdppc ~  l_subsidy_rgdppc  +  lag_gtax_pc_cons    + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party  +  l_sl_neigh_greentax   +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub #  + ccode 
labels<-list("tax2", "subsidy2") 
system<-list(tax2, subsidy2)                                                       
inst1<- ~  lag_gtax_pc_cons+ rgdppc + rgdppcsq + unemp+ unempsq + actual_economic_flows  +  rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax    + ccode 
inst2<- ~   l_subsidy_rgdppc + rgdppc + rgdppcsq + unemp+ unempsq +  actual_economic_flows   + rile_const_w +  env_const_w +  env_const_w_sq + green_party + en_prod_rgdp + sl_gled_sub +  sl_igo_sub + sl_dtrade_sub +    ccode                                                                                    
inst<-list(inst1, inst2)
fit3sls_closecountries_a <-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_closecountries_a)  

#  2SLS, simultanous equations - conditional osmosis of environmental policies (by acual economic flows)
tax3   <- greentax_pc_cons  ~ lag_gtax_pc_cons + l_subsidy_rgdppc +  i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party +  l_sl_neigh_subsidy_rgdppc +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  + I( l_sl_neigh_subsidy_rgdppc*actual_economic_flows)    # 
subsidy3  <- subsidy_rgdppc ~  l_subsidy_rgdppc   + lag_gtax_pc_cons   + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party  +  l_sl_neigh_greentax + sl_gled_sub +  sl_igo_sub + sl_dtrade_sub  + I( l_sl_neigh_greentax*actual_economic_flows)  #  
labels<-list("tax3", "subsidy3") 
system<-list(tax3, subsidy3)                                                       
inst1<- ~ lag_gtax_pc_cons+ rgdppc + rgdppcsq+ unemp+ unempsq     + en_prod_rgdp  +   sl_gled_tax +  sl_igo_tax + sl_dtrade_tax+     ccode 
inst2<- ~  l_subsidy_rgdppc + lag_gtax_pc_cons + rgdppc + rgdppcsq+ unemp+ unempsq  + en_prod_rgdp +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub  + ccode                                                                                   
inst<-list(inst1, inst2)
fit3sls_closecountries_b <-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_closecountries_b)  

##############
# Figure A.2 #
##############

tax.fig <- lm(greentax_pc_cons  ~  lag_gtax_pc_cons  + l_subsidy_rgdppc   + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp  + l_sl_neigh_subsidy_rgdppc +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax + I( l_sl_neigh_subsidy_rgdppc*actual_economic_flows)    , data=scr)
summary(tax.fig)
#visreg2d(tax.fig, y="l_sl_neigh_subsidygdp_st", x="actual_economic_flows",plot.type="image", main="conditional effects on green tax levels", ylab="subsidy spatial lag", xlab="actual economic flows", cex.lab=.9, cex.axis=.5, cex.main=.9, cex.sub=.6)
visreg2d(tax.fig, y="l_sl_neigh_subsidy_rgdppc", x="actual_economic_flows",plot.type="persp", zlab="green tax", ylab="subsidy spatial lag", xlab="actual economic flows", cex.lab=.9, cex.axis=.8, cex.main=.5, cex.sub=.6)

###################################
# All distances: Read the data  #
###################################

# Read the data:
scr <- read.dta("~/Dropbox/.../gkm_replication_isq/green taxes and subsidies/OECD subsidies.dta")
scr$ccode<-as.factor(as.character(scr$ccode))
scr$year<-as.factor(as.character(scr$year))

######################################################
#  All distances, OLS: first 3 columns of Table A.5  #
######################################################

tax1a     <- greentax_pc_cons  ~   lag_gtax_pc_cons + l_subsidy_rgdppc   + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp  +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  +   ccode    
subsidy1a     <- subsidy_rgdppc ~   l_subsidy_rgdppc + lag_gtax_pc_cons    + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp   +     sl_gled_sub +  sl_igo_sub + sl_dtrade_sub  +    ccode  
labels<-list("tax1a", "subsidy1a") 
system<-list(tax1a, subsidy1a) 
fitolsa <- systemfit( system, data=scr )
summary(fitolsa)
out <- capture.output(summary(fitolsa))

tax1b    <- greentax_pc_cons  ~   lag_gtax_pc_cons + l_subsidy_rgdppc    + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp  +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  + l_sl_gled_sub    + ccode
subsidy1b    <- subsidy_rgdppc ~   l_subsidy_rgdppc + lag_gtax_pc_cons    + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp   +     sl_gled_sub +  sl_igo_sub + sl_dtrade_sub + l_sl_gled_tax     +ccode
labels<-list("tax1b", "subsidy1b") 
system<-list(tax1b, subsidy1b) 
fitolsb <- systemfit( system, data=scr )
summary(fitolsb)
out <- capture.output(summary(fitolsb))

tax1c   <- greentax_pc_cons  ~   lag_gtax_pc_cons + l_subsidy_rgdppc    + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp  +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  + l_sl_gled_sub +    I( l_sl_gled_sub*actual_economic_flows)     + ccode
subsidy1c    <- subsidy_rgdppc ~   l_subsidy_rgdppc + lag_gtax_pc_cons    + rgdppc + rgdppcsq + unemp + unempsq + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp   +     sl_gled_sub +  sl_igo_sub + sl_dtrade_sub + l_sl_gled_tax  + I( l_sl_gled_tax*actual_economic_flows)    +ccode
labels<-list("tax1c", "subsidy1c") 
system<-list(tax1c, subsidy1c) 
fitolsc <- systemfit( system, data=scr )
summary(fitolsc)
out <- capture.output(summary(fitolsc))

###############
# Figure A.3  #
###############

tax_simple_gdp     <- lm(greentax_pc_cons  ~   lag_gtax_pc_cons + l_subsidy_rgdppc + rgdppc + rgdppcsq + unemp + unempsq  +  i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp  +  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax + l_sl_gled_sub  + I( l_sl_gled_sub*actual_economic_flows)  , data =scr)
#visreg2d(tax_simple_gdp, y="l_sl_gled_sub", x="actual_economic_flows",plot.type="image", main="conditional effects on green tax levels", ylab="subsidy spatial lag", xlab="actual economic flows", cex.lab=.9, cex.axis=.5, cex.main=.96, cex.sub=.75)
visreg2d(tax_simple_gdp, y="l_sl_gled_sub", x="actual_economic_flows",plot.type="persp", zlab="green tax", ylab="subsidy spatial lag", xlab="actual economic flows", cex.lab=.9, cex.axis=.8, cex.main=.5, cex.sub=.6)

####################################################
# Additional table: all distances, 3SLS w/ GDP p.c. #
####################################################

tax2a   <- greentax_pc_cons  ~ lag_gtax_pc_cons +   l_subsidy_rgdppc  + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party  + sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  # 
subsidy2a  <- subsidy_rgdppc ~ l_subsidy_rgdppc  +  lag_gtax_pc_cons   + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party     +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub# 
labels<-list("tax2a", "subsidy2a") 
system<-list(tax2a, subsidy2a)                                                       
inst1<- ~   rgdppc + rgdppcsq + unemp+ unempsq + actual_economic_flows +    i_itax_per_cap +  rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp + sl_gled_tax +  sl_igo_tax + sl_dtrade_tax    + ccode 
inst2<- ~   rgdppc + rgdppcsq + unemp+ unempsq +  actual_economic_flows +   i_itax_per_cap  + rile_const_w +  env_const_w +  env_const_w_sq + green_party + en_prod_rgdp +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub+    ccode                                                                                    
inst<-list(inst1, inst2)
fit3sls_a<-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_a)  

tax2b   <- greentax_pc_cons  ~ lag_gtax_pc_cons +   l_subsidy_rgdppc  + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party  + l_sl_gled_sub + sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  # 
subsidy2b  <- subsidy_rgdppc ~ l_subsidy_rgdppc  +  lag_gtax_pc_cons   + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party  +  l_sl_gled_tax   +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub# 
labels<-list("tax2b", "subsidy2b") 
system<-list(tax2b, subsidy2b)                                                       
inst1<- ~   rgdppc + rgdppcsq + unemp+ unempsq + actual_economic_flows +    i_itax_per_cap +  rile_const_w + env_const_w +  env_const_w_sq + green_party + en_prod_rgdp + sl_gled_tax +  sl_igo_tax + sl_dtrade_tax    + ccode 
inst2<- ~   rgdppc + rgdppcsq + unemp+ unempsq +  actual_economic_flows +   i_itax_per_cap  + rile_const_w +  env_const_w +  env_const_w_sq + green_party + en_prod_rgdp +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub+    ccode                                                                                    
inst<-list(inst1, inst2)
fit3sls_b<-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_b)  

#  3SLS, simultanous equations - (by acual economic flows)

tax2c  <- greentax_pc_cons  ~ lag_gtax_pc_cons + l_subsidy_rgdppc +  i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party +  l_sl_gled_sub + sl_gled_tax +  sl_igo_tax + sl_dtrade_tax  + I( l_sl_gled_sub*actual_economic_flows)    # 
subsidy2c  <- subsidy_rgdppc ~  l_subsidy_rgdppc   + lag_gtax_pc_cons  + i_itax_per_cap + actual_economic_flows + rile_const_w + env_const_w +  env_const_w_sq + green_party  +  l_sl_gled_tax +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub + I( l_sl_gled_tax*actual_economic_flows)  # 
labels<-list("tax2c", "subsidy2c") 
system<-list(tax2c, subsidy2c)                                                       
inst1<- ~  rgdppc + rgdppcsq+ unemp+ unempsq  +  actual_economic_flows +   i_itax_per_cap  + rile_const_w +  env_const_w +  env_const_w_sq + green_party + en_prod_rgdp+  sl_gled_tax +  sl_igo_tax + sl_dtrade_tax+     ccode 
inst2<- ~   rgdppc + rgdppcsq+ unemp+ unempsq +  actual_economic_flows +   i_itax_per_cap  + rile_const_w +  env_const_w +  env_const_w_sq + green_party + en_prod_rgdp +  sl_gled_sub +  sl_igo_sub + sl_dtrade_sub + ccode                                                                                   
inst<-list(inst1, inst2)
fit3sls_c<-systemfit( system,"3SLS",  inst=inst,  data = scr )  
summary(fit3sls_c)  
